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Abstract. Interstellar is the first Hollywood movie to attempt depicting a black hole 
as it would actually be seen by somebody nearby. For this, our team at Double Negative 
Visual Effects, in collaboration with physicist Kip Thorne, developed a code called 
DNGR (Double Negative Gravitational Renderer) to solve the equations for ray-bundle 
(light-beam) propagation through the curved spacetime of a spinning (Kerr) black hole, 
and to render IMAX-quality, rapidly changing images. Our ray-bundle techniques were 
crucial for achieving IMAX-quality smoothness without flickering; and they differ from 
physicists’ image-generation techniques (which generally rely on individual light rays 
rather than ray bundles), and also differ from techniques previously used in the film 
industry’s CGI community. 

This paper has four purposes: (i) To describe DNGR for physicists and CGI 
practitioners, who may find interesting and useful some of our unconventional 
techniques, (ii) To present the equations we use, when the camera is in arbitrary 
motion at an arbitrary location near a Kerr black hole, for mapping light sources to 
camera images via elliptical ray bundles, (iii) To describe new insights, from DNGR, 
into gravitational lensing when the camera is near the spinning black hole, rather than 
far away as in almost all prior studies; we focus on the shapes, sizes and influence 
of caustics and critical curves, the creation and annihilation of stellar images, the 
pattern of multiple images, and the influence of almost-trapped light rays, and we 
find similar results to the more familiar case of a camera far from the hole, (iv) To 
describe how the images of the black hole Gargantua and its accretion disk, in the 
movie Interstellar, were generated with DNGR—including, especially, the influences 
of (a) colour changes due to doppler and gravitational frequency shifts, (b) intensity 
changes due to the frequency shifts, (c) simulated camera lens flare, and (d) decisions 
that the film makers made about these influences and about the Gargantua’s spin, with 
the goal of producing images understandable for a mass audience. There are no new 
astrophysical insights in this accretion-disk section of the paper, but disk novices may 
find it pedagogically interesting, and movie buffs may find its discussions of Interstellar 
interesting. 
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1. Introduction and overview 
1.1. Previous research and visualizations 

At a summer school in Les Houches France in summer 1972, James Bardeen |T], building 
on earlier work of Brandon Carter [2], initiated research on gravitational lensing by 
spinning black holes. Bardeen gave a thorough analytical analysis of null geodesics 
(light-ray propagation) around a spinning black hole; and, as part of his analysis, he 
computed how a black hole’s spin affects the shape of the shadow that the hole casts 
on light from a distant star held. The shadow bulges out on the side of the hole moving 
away from the observer, and squeezes inward and hattens on the side moving toward 
the observer. The result, for a maximally spinning hole viewed from afar, is a D-shaped 
shadow; cf. Figure |^below. (When viewed up close, the shadow’s hat edge has a shallow 
notch cut out of it, as hinted by Figure below.) 

Despite this early work, gravitational lensing by black holes remained a backwater 
of physics research until decades later, when the prospect for actual observations brought 
it to the fore. 

There were, we think, two especially memorable accomplishments in the backwater 
era. The hrst was a 1978 simulation of what a camera sees as it orbits a non-spinning 
black hole, with a star held in the background. This simulation was carried out by 
Leigh Palmer, Maurice Pryce and Bill Unruh [3] on an Evans and Sutherland Vector 
graphics display at Simon Fraser University. Palmer, Pryce and Unruh did not publish 
their simulation, but they showed a him clip from it in a number of lectures in that era. 
The nicest modern-era him clip of this same sort that we know of is by Alain Riazuelo 
(contained in his DVD |3] and available on the web at m ; see Figure]^ and associated 
discussion below. And see [6] for an online application by Thomas Muller and Daniel 
Weiskopf for generating similar him clips. Also of much interest in our modern era are 
him clips by Andrew Hamilton [7] of what a camera sees when falling into a nonspinning 
black hole; these have been shown at many planetariums, and elsewhere. 

The other most memorable backwater-era accomplishment was a black and white 
simulation by Jean-Pierre Luminet [8] of what a thin accretion disk, gravitationally 
lensed by a nonspinning black hole, would look like as seen from far away but close 
enough to resolve the image. In Figure [T^ below, we show a modern-era colour version 
of this, with the camera close to a fast-spinning black hole. 

Gravitational lensing by black holes began to be observationally important in the 
1990s. 

Kevin Rauch and Roger Blandford jH] recognised that, when a hot spot, in a 
black hole’s accretion disk or jet, passes through caustics of the Earth’s past light cone 
(caustics produced by the hole’s spacetime curvature), the brightness of the hot spot’s 
X-rays will undergo sharp oscillations with informative shapes. This has motivated a 
number of quantitative studies of the Kerr metric’s caustics; see, especially Bdomi] 
and references therein. 

These papers’ caustics are relevant for a source near the black hole and an observer 
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far away, on Earth—in effect, on the black hole’s “celestial sphere” at radius r = cx). 
In our paper, by contrast, we are interested in light sources that are usually on the 
celestial sphere and an observer or camera near the black hole. For this reversed case, 
we shall discuss the relevant caustics in Sections 13.31 and 13.41 This case has been little 
studied, as it is of primarily cultural interest (“everyone” wants to know what it would 
look like to live near a black hole, but nobody expects to make such observations in 
his or her lifetime), and of science-hction interest. Our paper initiates the detailed 
study of this culturally interesting case; but we leave a full, systematic study to future 
research. Most importantly, we keep our camera outside the ergosphere — by contrast 
with Alain Riazuelo’s [ 12 ] and Hamilton’s [ 7 ] recent simulations with cameras deep 
inside the ergosphere and even plunging through the horizon||] 

In the 1990s astrophysicists began to envision an era in which very long baseline 
interferometry would make possible the imaging of black holes — specihcally, their 
shadows and their accretion disks. This motivated visualizations, with ever increasing 
sophistication, of accretion disks around black holes: modern variants of Luminet’s [8] 
pioneering work. See, especially, Fukue and Yokoyama [13], who added colours to the 
disk; Viergutz [TT], who made his black hole spin, treated thick disks, and produced 
particularly nice and interesting coloured images and included the disk’s secondary 
image which wraps under the black hole; Marck [IS], who laid the foundations for a 
lovely movie now available on the web [I6] with the camera moving around close to the 
disk, and who also included higher-order images, as did Fanton et. al. im and Beckwith 
and Done [IH]. See also papers cited in these articles. 

In the 2000s astrophysicists have focused on perfecting the mm-interferometer 
imaging of black-hole shadows and disks, particularly the black hole at the centre of our 
own Milky Way galaxy (Sgr A*). See, e.g., the 2000 feasibility study by Falcke, Melia 
and Agol [12]. See also references on the development and exploitation of GRMHD 
(general relativistic magnetohydrodynamical) simulation codes for modelling accretion 
disks like that in Sgr A* [20l |2T1[22]; and references on detailed GRMHD models of Sgr 
A* and the models’ comparison with observations [2HII241I25112B] . This is culminating in 
a mm interferometric system called the Event Horizon Telescope 123, which is beginning 
to yield interesting observational results though not yet full images of the shadow and 
disk in Sgr A*. 

All the astrophysical visualizations of gravitational lensing and accretion disks 
described above, and all others that we are aware of, are based on tracing huge numbers 
of light rays through curved spacetime. A primary goal of today’s state-of-the-art, 
astrophysical ray-tracing codes (e.g., the Ghan, Psaltis and Ozel’s massively parallel, 
GPU-based code GRay [28]) is very fast throughput, measured, e.g., in integration steps 
per second; the spatial smoothness of images has been only a secondary concern. For our 
Interstellar work, by contrast, a primary goal is smoothness of the images, so flickering 
is minimised when objects move rapidly across an IMAX screen; fast throughput has 

f We also have done such simulations but not in enough detail to reveal much new. 
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been only a secondary concern. 

With these different primary goals, in our own code, called DNGR, we have been 
driven to employ a different set of visualization techniques from those of the astrophysics 
community—techniques based on propagation of ray bundles (light beams) instead of 
discrete light rays, and on carefully designed spatial filtering to smooth the overlaps of 
neighbouring beams; see Section]^ and Appendix A Although, at Double Negative, we 
have many GPU-based workstations, the bulk of our computational work is done on a 
large compute cluster (the Double Negative render-farm) that does not employ GPUs. 

In [Appendix A.7| we shall give additional comparisons of our DNGR code with 
astrophysical codes and with other film-industry GGI codes. 


1.2. This paper 

Our work on gravitational lensing by black holes began in May 2013, when Ghristopher 
Nolan asked us to collaborate on building realistic images of a spinning black hole and 
its disk, with IMAX resolution, for his science fiction movie Interstellar. We saw this not 
only as an opportunity to bring realistic black holes into the Hollywood arena, but also 
an opportunity to create a simulation code capable of exploring a black hole’s lensing 
with a level of image smoothness and dynamics not previously available. 

To achieve IMAX quality (with 23 million pixels per image and adequately 
smooth transitions between pixels), our code needed to integrate not only rays (photon 
trajectories) from the light source to the simulated camera, but also bundles of rays 
(light beams) with filtering to smooth the beams’ overlap; see Section]^, Appendix A.2 


and [Appendix A.3[ And because the camera would sometimes be moving with speeds 


that are a substantial fraction of the speed of light, our code needed to incorporate 
relativistic aberration as well as Doppler shifts and gravitational redshifts. 

Thorne, having had a bit of experience with this kind of stuff, put together a step- 
by-step prescription for how to map a light ray and ray bundle from the light source (the 


celestial sphere or an accretion disk) to the camera’s local sky; see Appendix A.l and 


Appendix A.2 He implemented his prescription in Mathematica to be sure it produced 
images in accord with others’ prior simulations and his own intuition. He then turned 
his prescription over to our Double Negative team, who created the fast, high-resolution 
code DNGR that we describe in Section and Appendix A and created the images to 
be lensed: fields of stars and in some cases also dust clouds, nebulae, and the accretion 
disk around Interstellar’s black hole, Gargantua. 

In Section we use this code to explore a few interesting aspects of the lensing 
of a star field as viewed by a camera that moves around a fast-spinning black hole 
(0.999 of maximal—a value beyond which natural spin-down processes become strong 
m)- We compute the intersection of the celestial sphere (r = oo) with the primary, 
secondary, and tertiary caustics of the camera’s past light cone and explore how the 
sizes and shapes of the resulting caustic curves change as the camera moves closer to 
the hole and as we go from the primary caustic to higher-order caustics. We compute 
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the images of the hrst three caustics on the camera’s local sky (the hrst three critical 
curves), and we explore in detail the creation and annihilation of stellar-image pairs, on 
the secondary critical curve, when a star passes twice through the secondary caustic. 
We examine the tracks of stellar images on the camera’s local sky as the camera orbits 
the black hole in its equatorial plane. And we discover and explain why, as seen by a 
camera orbiting a fast spinning black hole in our galaxy, there is just one image of the 
galactic plane between adjacent critical curves near the poles of the black hole’s shadow, 
but there are multiple images between critical curves near the flat, equatorial shadow 
edge. The key to this from one viewpoint is the fact that higher-order caustics wrap 
around the celestial sphere multiple times, and from another viewpoint the key is light 
rays temporarily trapped in prograde, almost circular orbits around the black hole. By 
placing a checkerboard of paint swatches on the celestial sphere, we explore in detail the 
overall gravitational lensing patterns seen by a camera near a fast-spinning black hole 
and the influence of aberration due to the camera’s motion. 

Whereas most of Section on lensing of stellar images is new, in the context of a 
camera near the hole and stars far away, our Section on lensing of an accretion disk, 
retreads largely known ground. But it does so in the context of the movie Interstellar, 
to help readers understand the movie’s black-hole images, and does so in a manner that 
may be pedagogically interesting. We begin with a picture of a gravitationally lensed 
disk made of equally spaced paint swatches. This picture is useful for understanding 
the multiple images of the disk that wrap over and under and in front of the hole’s 
shadow. We then replace the paint-swatch disk by a fairly realistic and rather thin disk 
(though one constructed by Double Negative artists instead of by solving astrophysicists’ 
equations for thin accretion disks [ 32 ] )• We hrst compute the lensing of our semi- 
realistic disk but ignore Doppler shifts and gravitational redshifts, which we then turn 
on pedagogically in two steps: the colour shifts and then the intensity shifts. We discuss 
why Christopher Nolan and Paul Franklin chose to omit the Doppler shifts in the movie, 
and chose to slow the black hole’s spin from that, a/M ~ 1, required to explain the huge 
time losses in Interstellar, to a/M = 0.6. And we then discuss and add simulated lens 
hare (light scattering and dihraction in the camera’s lenses) for an IMAX camera that 
observes the disk - something an astrophysicist would not want to do because it hides 
the physics of the disk and the lensed galaxy beyond it, but that is standard in movies, 
so computer generated images will have continuity with images shot by real cameras. 

Finally, in Section]^ we summarise and we point to our use of DNGR to produce 
images of gravitational lensing by wormholes. 

Throughout we use geometrized units in which G (Newton’s gravitation constant) 
and c (the speed of light) are set to unity, and we use the MTW sign conventions [33] . 

2. DNGR: our Double Negative Gravitational Renderer code 

Our computer code for making images of what a camera would see in the vicinity of a 
black hole or wormhole is called the Double Negative Gravitational Renderer, or DNGR 
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— which obviously can also be interpreted as the Double Negative General Relativistic 
code. 


2.1. Ray Tracing 


The ray tracing part of DNGR produces a map from the celestial sphere (or the surface 
of an accretion disk) to the camera’s local sky. More specifically (see 
details and Figurefor the ray-tracing geometry): 


Appendix A.l for 


(i) In DNGR, we adopt Boyer-Lindquist coordinates (t, r, 6*, 0) for the black hole’s Kerr 
spacetime. At each event we introduce the locally non-rotating observer, also called 
the fiducial observer or FIDO in the Membrane Paradigm [31]: the observer whose 
4-velocity is orthogonal to the surfaces of constant t, the Kerr metric’s space slices. 
We regard the FIDO as at rest in space, and give the FIDO orthonormal basis 
vectors, {er,eg,e^}, that point along the spatial coordinate lines. 



Figure 1: The mapping of the camera’s local sky (^c*, 0cs) onto the celestial sphere 
{O', (/>') via a backward directed light ray; and the evolution of a ray bundle, that is 
circular at the camera, backward along the ray to its origin, an ellipse on the celestial 
sphere. 







Gravitational Lensing by Black Holes in Astrophysics and in Interstellar 


7 


(ii) We specify the camera’s coordinate location {rc,9c,(j)c)] its direction of motion 
relative to the FIDO there, a unit vector B in the camera’s reference frame; and 
the camera’s speed f3 relative to the FIDO. 

(hi) In the camera’s reference frame, we set up a right-handed set of three orthonormal 
basis vectors {ex,ey,ez}, with = B along the direction of the camera’s motion, 
e^, perpendicular to ey and in the plane spanned by and e^, and orthogonal to 
ea, and By. See Figure And we then set up a spherical polar coordinate system 
{0cs,4>cs} for the camera’s local sky (i.e. for the directions of incoming light rays) 
in the usual manner dictated by the camera’s Cartesian basis vectors {ex,ey,ez}. 

(iv) For a ray that originates on the celestial sphere (at r = oo), we denote the Boyer- 
Lindquist angular location at which it originates by 

(v) We integrate the null geodesic equation to propagate the ray from the camera to 
the celestial sphere, thereby obtaining the map {9'{9cs,(f)cs),4>'{9cs,((>cs)} of points 
on the camera’s local sky to points on the celestial sphere. 


(vi) 


(vii) 


If the ray originates on the surface of an accretion disk, we integrate the null 
geodesic equation backward from the camera until it hits the disk’s surface, and 
thereby deduce the map from a point on the disk’s surface to one on the camera’s 


sky. For more details on this case, see Appendix A.6 


We also compute, using the relevant Doppler shift and gravitational redshift, the 
net frequency shift from the ray’s source to the camera, and the corresponding net 
change in light intensity. 


2.2. Ray-Bundle (Light-Beam) Propagation 


DNGR achieves its IMAX-quality images by integrating a bundle of light rays (a light 
beam) backward along the null geodesic from the camera to the celestial sphere using a 
slightly modihed variant of a procedure formulated in the 1970s by Serge Pineault and 
Robert Roeder [3Sl EE] • This procedure is based on the equation of geodesic deviation 
and is equivalent to the optical scalar equations EH that have been widely used by 
astrophysicists in analytical (but not numerical) studies of gravitational lensing; see 
references in Section 2.3 of [3H]. Our procedure, in brief outline, is this (see Figure [^; 
for full details, see [Appendix A.2 


(i) In DNGR, we begin with an initially circular (or sometimes initially elliptical) 
bundle of rays, with very small opening angle, centred on a pixel on the camera’s 
sky. 

(ii) We integrate the equation of geodesic deviation backward in time along the bundle’s 
central ray to deduce the ellipse on the celestial sphere from which the ray bundle 
comes. More specihcally, we compute the angle /i that the ellipse’s major axis 
makes with the celestial sphere’s e^/ direction, and the ellipse’s major and minor 
angular diameters 5+ and on the celestial sphere. 
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(ill) We then add up the spectrum and intensity of all the light emitted from within that 
ellipse; and thence, using the frequency and intensity shifts that were computed by 
ray tracing, we deduce the spectrum and intensity of the light arriving in the chosen 
camera pixel. 


2.3. Filtering, Implementation, and Code Characteristics 


Novel types of filtering are key to generating our IMAX-quality images for movies. In 
DNGR we use spatial hltering to smooth the interfaces between beams (ray bundles), 
and temporal filtering to make dynamical images look like they were filmed with a movie 
camera. For details, see [Appendix A.3 


In Appendix A.4 we describe some details of our DNGR implementation of the 
ray-tracing, ray-bundle, and hltering equations; in Appendix A.5 we describe some 
characteristics of our code and of Double Negative’s Linux-based render-farm on which 


we do our computations; in Appendix A.6 we describe our DNGR modelling of accretion 
disks; and in Appendix A.7| we briehy compare DNGR with other hlm-industry GGI 
codes and state-of-the-art astrophysical simulation codes. 


3. Lensing of a star field as seen by a moving camera near a black hole 

3.1. Nonspinning black hole 

In this subsection we review well known features of gravitational lensing by a 
nonspinning (Schwarzschild) black hole, in preparation for discussing the same things 
for a fast-spinning hole. 

We begin, pedagogically, with a still from a him clip by Alain Riazuelo [5], Figure 
The camera, at radius Vc = SOM (where M is the black hole’s mass) is moving in a 
circular geodesic orbit around the black hole, with a star held on the celestial sphere. 
We focus on two stars, which each produce two images on the camera’s sky. We place 
red circles around the images of the brighter star and yellow diamonds around those of 
the dimmer star. As the camera orbits the hole, the images move around the camera’s 
sky along the red and yellow curves. 

Images outside the Einstein ring (the violet circle) move rightward and dehect 
away from the ring. These are called primary images. Images inside the Einstein ring 
{secondary images) appear, in the him clip, to emerge from the edge of the black hole’s 
shadow, loop leftward around the hole, and descend back into the shadow. However, 
closer inspection with higher resolution reveals that their tracks actually close up along 
the shadow’s edge as shown in the hgure; the close-up is not seen in the him clip because 
the images are so very dim along the inner leg of their tracks. At all times, each star’s 
two images are on opposite sides of the shadow’s centre. 

This behaviour is generic. Every star (if idealised as a point source of light), except 
a set of measure zero, has two images that behave in the same manner as the red and 
yellow ones. Outside the Einstein ring, the entire primary star held hows rightward. 
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Figure 2: Gravitational lensing of a star field by a nonspinning black hole, as seen by a 
camera in a circular geodesic orbit at radius r = SOM. Picture courtesy Alain Riazuelo, 
from his him clip [5]; coloured markings by us. 


dehecting around the ring; inside the ring, the entire secondary star held loops leftward, 
conhned by the ring then back rightward along the shadow’s edge. (There actually are 
more, unseen, images of the star held, even closer to the shadow’s edge, that we shall 
discuss in Section 3.2[ ) 

As is well known, this behaviour is easily understood by tracing light rays from the 
camera to the celestial sphere; see Figure 

The Einstein ring is the image, on the camera’s sky, of a point source that is on 
the celestial sphere, diametrically opposite the camera; i.e., at the location indicated 
by the red dot and labeled “Caustic” in Figure Light rays from that caustic point 
generate the purple ray surface that converges on the camera, and the Einstein ring is 
the intersection of that ray surface with the camera’s local sky. 

[The caustic point (red dot) is actually the intersection of the celestial sphere with 
a caustic line (a one-dimensional sharp edge) on the camera’s past light cone. This 
caustic line extends radially from the black hole’s horizon to the caustic point.] 

The hgure shows a single star (black dot) on the celestial sphere and two light rays 
that travel from that star to the camera, gravitationally deflecting around opposite sides 
of the black hole. One of these rays, the primary one, arrives at the camera outside the 
Einstein ring; the other, secondary ray, arrives inside the Einstein ring. 

Because the caustic point and the star on the celestial sphere both have dimension 
zero, as the camera moves, causing the caustic point to move relative to the star, there 
is zero probability for it to pass through the star. Therefore, the star’s two images will 
never cross the Einstein ring; one will remain forever outside it and the other inside—and 
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Ray Surface From Caustic, through Einstein Ring to Camera 



Figure 3: Light rays around a Schwarzschild Black hole: geometric construction for 
explaining Figure 


similarly for all other stars in the star held. 

However, if a star with hnite size passes close to the ring, the gravitational lensing 
will momentarily stretch its two images into lenticular shapes that hug the Einstein ring 
and will produce a great, temporary increase in each image’s energy hux at the camera 
due to the temporary increase in the total solid angle subtended by each lenticular 
image. This increase in hux still occurs when the star’s actual size is too small for its 
images to be resolved, and also in the limit of a point star. For examples, see Riazuelo’s 
him clip [S]. 

(Large amplihcations of extended images are actually seen in Nature, for example in 
the gravitational lensing of distant galaxies by more nearby galaxies or galaxy clusters; 
see, e.g., [39].) 

3.2. Fast spinning black hole: introduction 

For a camera orbiting a spinning black hole and a star held (plus sometimes dust clouds 
and nebulae) on the celestial sphere, we have carried out a number of simulations with 
our code DNGR. We show a few him clips from these simulations at |10|. Figure]^ is 
a still from one of those him clips, in which the hole has spin a/M = 0.999 (where a is 
the hole’s spin angular momentum per unit mass and M is its mass), and the camera 
moves along a circular, equatorial, geodesic orbit at radius Vc = 6.03M. 

In the hgure we show in violet two critical curves—analogs of the Einstein ring for 
a nonspinning black hole. These are images, on the camera sky, of two caustic curves 






Gravitational Lensing by Black Holes in Astrophysics and in Interstellar 


11 



Figure 4: Gravitational lensing of a star field by a black hole with spin parameter 
a/M = 0.999, as seen by a camera in a circular, equatorial geodesic orbit at radius 
Vc = 6.03M. The red curves are the trajectories of primary images, on the camera’s sky, 
for stars at celestial-sphere latitudes 9' = O.IStt, 0.267r, 0.407r, O.OOtt, 0.747r, O.STtt. The 
yellow curves are the trajectories of secondary images for stars at 6' = O.SStt, O.Ohvr. The 
picture in this hgure is a still from our first film clip archived at [10] and is copyright 
©2015 Warner Bros. Entertainment Inc. Interstellar and all related characters and 
elements are trademarks of and ©Warner Bros. Entertainment Inc. (sl5). The full 
hgure appears in the second and later printings of The Science of Interstellar BH, and 
is used by permission of W. W. Norton & Company, Inc. This image may be used 
under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs 3.0 
(CC BY-NC-ND 3.0) license. Any further distribution of these images must maintain 
attribution to the author(s) and the title of the work, journal citation and DOI. You 
may not use the images for commercial purposes and if you remix, transform or build 
upon the images, you may not distribute the modihed images. 


that reside on the celestial sphere; see discussion below. 

We shall discuss in turn the region outside the secondary (inner) critical curve, and 
then the region inside. 

3.3. Fast-Spinning Hole: Outer region—outside the secondary critical curve 

As the camera moves through one full orbit around the hole, the stellar images in the 
outer region make one full circuit along the red and yellow curves and other curves 
like them, largely avoiding the two critical curves of Figure |4]— particularly the outer 
(primary) one. 
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We denote by five-pointed symbols four images of two very special stars: stars that 
reside where the hole’s spin axis intersects the celestial sphere, 6' = 0 and tt (0° and 
180°). These are analogs of the Earth’s star Polaris. By symmetry, these pole-star 
images must remain fixed on the camera’s sky as the camera moves along its circular 
equatorial orbit. Outside the primary (outer) critical curve, all northern-hemisphere 
stellar images (images with < 90°) circulate clockwise around the lower red pole- 
star image, and southern-hemisphere stellar images, counterclockwise around the upper 
red pole-star image. Between the primary and secondary (inner) critical curves, the 
circulations are reversed, so at the primary critical curve there is a divergent shear in 
the image flow. 

[For a nonspinning black hole (Figure above) there are also two critical curves, 
with the stellar-image motions confined by them: the Einstein ring, and a circular inner 
critical curve very close to the black hole’s shadow, that prevents the inner star tracks 
from plunging into the shadow and deflects them around the shadow so they close up.] 

3.3.1. Primary and secondary critical curves and their caustics After seeing these 
stellar-image motions in our simulations, we explored the nature of the critical curves 
and caustics for a camera near a fast-spinning black hole, and their influence. Our 
exploration, conceptually, is a rather straightforward generalisation of ideas laid out by 
Rauch and Blandford [9] and by Bozza [10]. They studied a camera or observer on the 
celestial sphere and light sources orbiting a black hole; our case is the inverse: a camera 
orbiting the hole and light sources on the celestial sphere. 

Just as the Einstein ring, for a nonspinning black hole, is the image of a caustic 
point on the celestial sphere—the intersection of the celestial sphere with a caustic 
line on the camera’s past light cone—so the critical curves for our spinning black hole 
are also images of the intersection of the celestial sphere with light-cone caustics. But 
the spinning hole’s light-cone caustics generically are 2-dimensional surfaces (folds) in 
the three-dimensional light cone, so their intersections with the celestial sphere are 
one-dimensional: they are closed caustic curves in the celestial sphere, rather than 
caustic points. The hole’s rotation breaks spherical symmetry and converts non-generic 
caustic points into generic caustic curves. (For this reason, theorems about caustics in 
the Schwarzschild spacetime, which are rather easy to prove, are of minor importance 
compared to results about generic caustics in the Kerr spacetime.) 

We have computed the caustic curves, for a camera near a spinning black hole, 
by propagating ray bundles backward from a fine grid of points on the camera’s local 
sky, and searching for points on the celestial sphere where the arriving ray bundle’s 
minor angular diameter passes through zero. All such points he on a caustic, and 
the locations on the camera sky where their ray bundles originate lie on critical curves. 

For a camera in the equatorial plane at radius Vc = 6.03M, Figure shows the 
primary and secondary caustic curves. These are images, in the celestial sphere, of the 
primary and secondary critical curves shown in Figure]^ 

The primary caustic is a very small astroid (a four-sided figure whose sides are fold 
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Figure 5: The primary and secondary caustics on the celestial sphere, for the past light 
cone of a camera moving along a circular, equatorial, geodesic orbit at radius Vc = 6.03M 
around a black hole with spin parameter a = 0.999M. As the camera moves, in the 
camera’s reference frame a star at 6' = OTOStt travels along the dashed-line path. 


caustics and meet in cusps). It spans just 0.45° in both the 0' and 9' directions. The 
secondary caustic, by contrast, is a large astroid: it extends over 263° in 0' and 72° in 9'. 
All stars within 36° of the equator encounter it as the camera, at Tc = 6.03M, orbits the 
black hole. This is similar to the case of a source near the black hole and a camera on 
the celestial sphere (far from the hole, e.g. on Earth) [101 [H]. There, also, the primary 
caustic is small and the secondary large. In both cases the dragging of inertial frames 
stretches the secondary caustic out in the 0 direction. 

3.3.2. Image creations and annihilations on critical curves Because the spinning hole’s 
caustics have finite cross sections on the celestial sphere, by contrast with the point 
caustics of a nonspinning black hole, stars, generically, can cross through them; see, 
e.g., the dashed stellar path in FigureAs is well known from the elementary theory 
of fold caustics (see, e.g.. Section 7.5 of |12]), at each crossing two stellar images, on 
opposite sides of the caustic’s critical curve, merge and annihilate; or two are created. 
And at the moment of creation or annihilation, the images are very bright. 

Figure]^ (two stills from a film clip at |10]) is an example of this. As the star in 
Figure |^, at polar angle 9' = O.bOSvr, travels around the celestial sphere relative to 
the camera, a subset of its stellar images travels around the red track of Figure]^ just 
below the black hole’s shadow. (These are called the star’s “secondary images” because 
the light rays that bring them to the camera have the same number of poloidal turning 
points, one—and equatorial crossings, one—as the light rays that map the secondary 
caustic onto the secondary critical curve; similarly these images’ red track is called the 
star’s “secondary track”.) At the moment of the left still, the star has just barely 
crossed the secondary caustic at point A of Figure [^, and its two secondary stellar 
images, 2 (inside the secondary critical curve) and 0^ 3 (outside it) have just been 
created at the point half way between and 0^3, where their red secondary track 
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Figure 6: Two stills from a film clip at [lO]). In the left still, images 2 and 3 have just 
been created as their star passed through caustic point A of Figure]^. In the right still, 
images 1 and 2 are about to annihilate as their star passes through caustic point B. 


crosses the secondary critical curve (Figure]^). In the meantime, stellar image is 
traveling slowly, alone, clockwise, around the track, outside the critical curve. Between 
the left and right stills, image ^2 travels the track counter clockwise and images 1 and 
3, clockwise. Immediately after the right still, the star crosses the secondary caustic at 
point B in Figure]^, and the two images (outside the critical curve) and #2 (inside 
it) annihilate at the intersection of their track with the critical curve (Figure]^). As the 
star in Figure continues on around the celestial sphere from point B to A, the lone 
remaining image on the track, image ^3, continues onward, clockwise, until it reaches 
the location of the hrst still, and two new images are created at the location between 
#2 and #3 of the hrst still; and so forth. 

Figure puts these in a broader context. It shows the tracks of all of the images 
of the star in Figure [^. Each image is labeled by its order n: the number of poloidal 
6 turning points on the ray that travels to it from its star on the celestial sphere; or, 
equally well (for our choice of a camera on the black hole’s equator), the number of 
times that ray crosses the equator 6 = 7r/2. The order-0 track is called the primary 
track, and (with no ray-equator crossings) it lies on the same side of the equator as its 
star; order-1 is the secondary track, and it lies on the opposite side of the equator from 
its star; order-2 is the tertiary track, on the same side of the equator as its star; etc. 

The primary track (order 0) does not intersect the primary critical curve, so a single 
primary image travels around it as the camera orbits the black hole. The secondary 
track (order 1) is the one depicted red in Figure and discussed above. It crosses 
the secondary critical curve twice, so there is a single pair creation event and a single 
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Figure 7: The tracks on the camera sky traveled by the multiple images of a single star, 
as the camera travels once around a black hole with a/M = 0.999. The camera’s orbit 
is a circular, equatorial geodesic with radius r^ = 6.03M. (a) For a star at latitude 
9' = O.OOtt (18° above the equatorial plane; essentially the same star as in Figure]^), 
(b) For a star at 6' = O.hlvr (1.8° above the equator). The tracks are labeled by the 
order of their stellar images (the number of poloidal, 6, turning points on the ray that 
brings an image to the camera). 


annihilation event; at some times there is a single secondary image on the track, and 
at others there are three. It is not clear to us whether the red secondary track crosses 
the tertiary critical curve (not shown); but if it does, there will be no pair creations 
or annihilations at the crossing points, because the secondary track and the tertiary 
critical curve are generated by rays with different numbers of poloidal turning points, 
and so the critical curve is incapable of influencing images on the track. The extension 
to higher-order tracks and critical curves, all closer to the hole’s shadow, should be clear. 
This pattern is qualitatively the same as when the light source is near the black hole 
and the camera far away, but in the hole’s equatorial plane m- 

And for stars at other latitudes the story is also the same; only the shapes of the 
tracks are changed. Figure is an example. It shows the tracks for a star just 1.8° 
above the black hole’s equatorial plane, at 6' = O.hlvr. 

The film clips at [ID] exhibit these tracks and images all together, and show a 
plethora of image creations and annihilations. Exploring these clips can be fun and 
informative. 
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3.4- Fast-Spinning Hole: Inner region—inside the secondary critical curve 

The version of DNGR that we used for Interstellar showed a surprisingly complex, 
hngerprint-like structure of gravitationally lensed stars inside the secondary critical 
curve, along the left side of the shadow. 

We searched for errors that might be responsible for it, and hnding none, we thought 
it real. But Alain Riazuelo |13] saw nothing like it in his computed images. Making 
detailed comparisons with Riazuelo, we found a bug in DNGR. When we corrected the 
bug, the complex pattern went away, and we got excellent agreement with Riazuelo 
(when using the same coordinate system), and with images produced by Andy Bohn, 
Francois Hebert and William Throwe using their Gornell/Galtech SXS imaging code |11]. 
Since the SXS code is so very different from ours (it is designed to visualize colliding 
black holes), that agreement gives us high conhdence in the results reported below. 

Fortunately, the bug we found had no noticeable impact on the images in 
Interstellar. 

With our debugged code, the inner region, inside the secondary critical curve, 
appears to be a continuation of the pattern seen in the exterior region. There is a third 
critical curve within the second, and there are signs of higher-order critical curves, all 
nested inside each other. These are most visible near the flattened edge of the black 
hole’s shadow on the side where the horizon’s rotation is toward the camera (the left 
side in this paper’s hgures). The dragging of inertial frames moves the critical curves 
outward from the shadow’s flattened edge, enabling us to see things that otherwise could 
only be seen with a strong zoom-in. 

3.4- 1. Critical curves and caustics forvc = 2.60M The higher-order critical curves and 
the regions between them can also be made more visible by moving the camera closer 
to the black hole’s horizon. In Figure we have moved the camera in to Vc = 2.6M 
with a/M = 0.999, and we show three nested critical curves. The primary caustic 
(the celestial-sphere image of the outer, primary, critical curve) is a tiny astroid, as at 
Vc = 6.03M (Figure]^. The secondary and tertiary caustics are shown in Figure]^ 

By comparing these three figures with each other, we see that (i) as we move the 
camera closer to the horizon, the secondary caustics wrap further around the celestial 
sphere (this is reasonable since frame dragging is stronger nearer the hole), and (ii) at 
fixed camera radius, each order caustic wraps further around the celestial sphere than 
the lower order ones. More specifically, for Vc = 2.60M, the secondary caustic (Figure 
[^) is stretched out to one full circuit around the celestial sphere compared to 2/3 of a 
circuit when the camera is at Vc = 6.03M (Figure]^, and the tertiary caustic (Figure 
[^) is stretched out to more than six circuits around the celestial sphere! 

The mapping of points between each caustic and its critical curve is displayed in 
him clips at [iQ]. For the secondary caustic at = 2.60M, we show a few points in that 
mapping in Figure |^. The left side of the critical curve (the side where the horizon 
is moving toward the camera) maps into the long, stretched-out leftward sweep of the 
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caustic. The right side maps into the caustic’s two unstretched right scallops. The same 
is true of other caustics and their critical curves. 

3.4-2. Multiple images for rc = 2.60M Returning to the gravitationally lensed star- 
held image in Figure [^: notice the series of images of the galactic plane (fuzzy white 
curves). Above and below the black hole’s shadow there is just one galactic-plane image 
between the primary and secondary critical curves, and just one between the secondary 
and tertiary critical curves. This is what we expect from the example of a nonspinning 
black hole. However, near the fast-spinning hole’s left shadow edge, the pattern is very 
different: three galactic-plane images between the primary and secondary critical curves, 
and eight between the secondary and tertiary critical curves. 

These multiple galactic-plane images are caused by the large sizes of the caustics— 
particularly their wrapping around the celestial sphere—and the resulting ease with 
which stars cross them, producing multiple stellar images. An extension of an argument 
by Bozza [10] [paragraph preceding his Eq. (17)] makes this more precise. (This 
argument will be highly plausible but not fully rigorous because we have not developed 
a sufficiently complete understanding to make it rigorous.) 

Consider, for concreteness, a representative galactic-plane star that lies in the black 
hole’s equatorial plane, inside the primary caustic, near the caustic’s left cusp; and ask 



Figure 8: (a) Black-hole shadow and three critical curves for a camera traveling on a 
circular, geodesic, equatorial orbit at radius Vc = 2.60M in the equatorial plane of a 
black hole that has spin a/M = 0.999. (b) Blowup of the equatorial region near the 
shadow’s flat left edge. The imaged star held is adapted from the Tycho-2 catalogue 
[15] of the brightest 2.5 million stars seen from Earth, so it shows multiple images of 
the galactic plane. 
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^cs, 0 

Figure 9: (a) The secondary caustic (red) on the celestial sphere and secondary critical 
curve (green) on the camera’s sky, for the black hole and camera of Figure]^ Points on 
each curve that are ray-mapped images of each other are marked by letters a, b, c, d. 
(b) The tertiary caustic and tertiary critical curve. 


how many images the star produces on the camera’s sky and where they he. To answer 
this question, imagine moving the star upward in 9' at hxed (j)' until it is above all the 
celestial-sphere caustics; then move it gradually back downward to its original, equatorial 
location. 

When above the caustics, the star produces one image of each order: A primary 
image (no poloidal turning points) that presumably will remain outside the primary 
critical curve when the star returns to its equatorial location; a secondary image (one 
poloidal turning point) that presumably will be between the primary and secondary 
critical curves when the star returns; a tertiary image between the secondary and tertiary 
critical curves; etc. 

When the star moves downward through the upper left branch of the astroidal 
primary caustic, it creates two primary images, one on each side of the primary critical 
curve. When it moves downward through the upper left branch of the secondary caustic 
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(Figure]^), it creates two secondary images, one on each side of the secondary caustic. 
And when it moves downward through the six sky-wrapped upper left branches of the 
tertiary caustic (Figure [^), it creates twelve tertiary images, six on each side of the 
tertiary caustic. And because the upper left branches of all three caustics map onto the 
left sides of their corresponding critical curves, all the created images will wind up on 
the left sides of the critical curves and thence the left side of the black hole’s shadow. 
And by symmetry, with the camera and the star both in the equatorial plane, all these 
images will wind up in the equatorial plane. 

So now we can count. In the equatorial plane to the left of the primary critical 
curve, there are two images: one original primary image, and one caustic-created 
primary image. These are to the left of the region depicted in Figure [^. Between the 
primary and secondary critical curves there are three images: one original secondary 
image, one caustic-created primary, and one caustic-created secondary image. These are 
representative stellar images in the three galactic-plane images between the primary and 
secondary critical curves of Figure And between the secondary and tertiary critical 
curves there are eight stellar images: one original tertiary, one caustic-created secondary, 
and six caustic-created tertiary images. These are representative stellar images in the 
eight galactic-plane images between the secondary and tertiary critical curves of Figure 

e 

This argument is not fully rigorous because: (i) We have not proved that every 
caustic-branch crossing, from outside the astroid to inside, creates an image pair rather 
than annihilating a pair; this is very likely true, with annihilations occurring when a 
star moves out of the astroid. (ii) We have not proved that the original order-n images 
wind up at the claimed locations, between the order-n and order-(n-|- 1) critical curves. 
A more thorough study is needed to pin down these issues. 


3.4-3. Checkerboard to elucidate the multiple-image pattern Figure pi)] is designed to 
help readers explore this multiple-image phenomenon in greater detail. There we have 


placed, on the celestial sphere, a checkerboard of paint swatches (Figure 10 r), with 
dashed lines running along the constant-latitude spaces between paint swatches, i.e., 
along the celestial-sphere tracks of stars. In Figure [To| 3 we show the gravitationally 
lensed checkerboard on the camera’s entire sky; and in Figure [TOj^ we show a blowup of 
the camera-sky region near the left edge of the black hole’s shadow. We have labeled 
the critical curves ICC, 2CC and 3CC for primary, secondary, and tertiary. 

The multiple images of lines of constant celestial-sphere longitude show up clearly 
in the blow-up, between pairs of critical curves; and the hgure shows those lines being 
stretched vertically, enormously, in the vicinity of each critical curve. The dashed lines 
(star-image tracks) on the camera’s sky show the same kind of pattern as we saw in 
Figure 


3.4-4- Multiple images explained by light-ray trapping The multiple images near the 
left edge of the shadow can also be understood in terms of the light rays that bring the 




Gravitational Lensing by Black Holes in Astrophysics and in Interstellar 


20 



0° 180° 360° 


Figure 10: (a) A checkerboard pattern of paint swatches placed on the celestial sphere of 
a black hole with spin a/M = 0.999. As the camera, moves around a circular, equatorial, 
geodesic orbit at radius Vc = 2.60M, stars move along horizontal dashed lines relative 
to the camera, (b) This checkerboard pattern as seen gravitationally lensed on the 
camera’s sky. Stellar images move along the dashed curves. The primary and secondary 
critical curves are labeled “ICC” and “2CC”. (c) Blowup of the camera’s sky near the 
left edge of the hole’s shadow; “3CC” is the tertiary critical curve. 
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Figure 11: A trapped, unstable, prograde photon orbit just outside the horizon of a 
black hole with spin a = 0.999M. The orbit, which has constant Boyer-Lindquist 
coordinate radius r, is plotted on a sphere, treating its Boyler-Lindquist coordinates 
as though they were spherical polar coordinates. 


stellar images to the camera. Those light rays travel from the celestial sphere inward to 
near the black hole, where they get temporarily trapped, for a few round-trips, on near 
circular orbits (orbits with nearly constant Boyer-Lindquist radius r), and then escape 
to the camera. Each such nearly trapped ray is very close to a truly (but unstably) 
trapped, constant-r ray such as that shown in Figure [II| These trapped rays (discussed 
in |16] and in Chapters 6 and 8 of [H]) wind up and down spherical strips with very 
shallow pitch angles. 

As the camera makes each additional prograde trip around the black hole, the image 
carried by each temporarily trapped mapping ray gets wound around the constant-r 
sphere one more time (i.e., gets stored there for one more circuit), and it comes out 
to the camera’s sky slightly closer to the shadow’s edge and slightly higher or lower in 
latitude. Correspondingly, as the camera moves, the star’s image gradually sinks closer 
to the hole’s shadow and gradually changes its latitude—actually moving away from the 
equator when approaching a critical curve and toward the equator when receding from 
a critical curve. This behaviour is seen clearly, near the shadow’s left edge, in the film 
clips at 


3.5. Aberration: Influence of the camera’s speed 


The gravitational lensing pattern is strongly influenced not only by the black hole’s spin 
and the camera’s location, but also by the camera’s orbital speed. We explore this in 
Figure [T^ where we show the gravitationally lensed paint-swatch checkerboard of Figure 


10 r for a black hole with spin a/M = 0.999, a camera in the equatorial plane at radius 


Tc = 2.60M, and three different camera velocities, all in the azimuthal direction: 
(a) camera moving along a prograde circular geodesic orbit [coordinate angular velocity 
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Figure 12: Influence of aberration, due to camera motion, on gravitational lensing by 
a black hole with a/M = 0.999. The celestial sphere is covered by the paint-swatch 
checkerboard of Figure [TO^ , the camera is at radius Vc = 2.60M and is moving in 
the azimuthal direction, and the camera speed is: (a) that of a prograde, geodesic, 
circular orbit (same as Figure [T^), (b) that of a zero-angular-momentum observer (a 
FIDO), and (c) at rest in the Boyer-Lindquist coordinate system. The coordinates are 
the same as in Figure 


10 D. 
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d(l)/dt = l/(a + r^/^) in the notation of Appendix A.l ; (b) camera moving along a zero- 
angular-momentum orbit [d(f)/dt = u, Eq. (A.2), which is speed —0.546c as measured by 
a circular, geodesic observer]; and (c) a static camera, i.e. at rest in the Boyer-Lindquist 
coordinate system [d(j)/dt = 0, which is speed —0.813c as measured by the circular, 
geodesic observer]. 

The huge differences in lensing pattern, for these three different camera velocities, 
are due, of course, to special relativistic aberration. (We thank Alain Riazuelo for 
pointing out to us that aberration effects should be huge.) For prograde geodesic 
motion (top picture), the hole’s shadow is relatively small and the sky around it, large. 
As seen in the geodesic reference frame, the zero-angular-momentum camera and the 
static camera are moving in the direction of the red/black dot—i.e., toward the right 
part of the external universe and away from the right part of the black-hole shadow— 
at about half and 4/5 the speed of light respectively. So the zero-angular-momentum 
camera (middle picture) sees the hole’s shadow much enlarged due to aberration, and 
the external universe shrunken; and the static camera sees the shadow enlarged so much 
that it encompasses somewhat more than half the sky (more than 27r steradians), and 
sees the external universe correspondingly shrunk. 

Despite these huge differences in lensing patterns, the multiplicity of images between 
critical curves is unchanged: still three images of some near-equator swatches between 
the primary and secondary critical curves, and eight between the secondary and tertiary 
critical curves. This is because the caustics in the camera’s past light cone depend only 
on the camera’s location and not on its velocity, so a point source’s caustic crossings are 
independent of camera velocity, and the image pair creations and annihilations along 
critical curves are independent of camera velocity. 


4. Lensing of an accretion disk 

4-1. Effects of lensing, colour shift, and brightness shift: a pedagogical discussion 

We have used our code, DNGR, to construct images of what a thin accretion disk 
in the equatorial plane of a fast-spinning black hole would look like, seen up close. 
For our own edification, we explored successively the influence of the bending of light 
rays (gravitational lensing), the influence of Doppler frequency shifts and gravitational 
frequency shifts on the disk’s colours, the influence of the frequency shifts on the 
brightness of the disk’s light, and the influence of lens flare due to light scattering and 
diffraction in the lenses of a simulated 65 mm IMAX camera. Although all these issues 
except lens flare have been explored previously, e.g. in [H [131 HU HU HH] and references 
therein, our images may be of pedagogical interest, so we show them here. We also 
show them as a foundation for discussing the choices that were made for Interstellar’s 
accretion disk. 
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Figure 13: Inset: Paint-swatch accretion disk with inner and outer radii r = 9.26M 
and r = 18.TOM before being placed around a black hole. Body: This paint-swatch 
disk, now in the equatorial plane around a black hole with a/M = 0.999, as viewed by 
a camera at r^ = 74.IM and 9c = 1.511 (86.56°), ignoring frequency shifts, associated 
colour and brightness changes, and lens flare. (Figure from The Science of Interstellar 
[H], used by permission of W. W. Norton & Company, Inc, and created by our Double 
Negative team, TM & (c)Warner Bros. Entertainment Inc. (sl5)). This image may be 
used under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs 
3.0 (CC BY-NC-ND 3.0) license. Any further distribution of these images must maintain 
attribution to the author(s) and the title of the work, journal citation and DOI. You 
may not use the images for commercial purposes and if you remix, transform or build 
upon the images, you may not distribute the modified images. 


4 . 1 . 1 . Gravitational lensing Figure 13 illustrates the influence of gravitational lensing 
(light-ray bending). To construct this image, we placed the infinitely thin disk shown 
in the upper left in the equatorial plane around a fast-spinning black hole, and we used 
DNGR to compute what the disk would look like to a camera near the hole and slightly 
above the disk plane. The disk consisted of paint swatches arranged in a simple pattern 
that facilitates identifying, visually, the mapping from the disk to its lensed images. 
We omitted frequency shifts and their associated colour and brightness changes, and 
also omitted camera lens flare; i.e., we (incorrectly) transported the light’s specific 
intensity A = dE/dtdAdQdu along each ray, unchanged, as would be appropriate in 
flat spacetime. Here E, t, A, D, and z/ are energy, time, area, solid angle, and frequency 
measured by an observer just above the disk or at the camera. 
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In the figure we see three images of the disk. The upper image swings around the 
front of the black hole’s shadow and then, instead of passing behind the shadow, it 
swings up over the shadow and back down to close on itself. This wrapping over the 
shadow has a simple physical origin: Light rays from the top face of the disk, which is 
actually behind the hole, pass up over the top of the hole and down to the camera due 
to gravitational light deflection; see Figure 9.8 of BH. This entire image comes from 
light rays emitted by the disk’s top face. By looking at the colours, lengths, and widths 
of the disk’s swatches and comparing with those in the inset, one can deduce, in each 
region of the disk, the details of the gravitational lensing. 

In Figure the lower disk image wraps under the black hole’s shadow and then 
swings inward, becoming very thin, then up over the shadow and back down and outward 
to close on itself. This entire image comes from light rays emitted by the disk’s bottom 
face: the wide bottom portion of the image, from rays that originate behind the hole, 
and travel under the hole and back upward to the camera; the narrow top portion, from 
rays that originate on the disk’s front underside and travel under the hole, upward on 
its back side, over its top, and down to the camera—making one full loop around the 
hole. 


There is a third disk image whose bottom portion is barely visible near the shadow’s 
edge. That third image consists of light emitted from the disk’s top face, that travels 
around the hole once for the visible bottom part of the image, and one and a half times 
for the unresolved top part of the image. 

In the remainder of this section we deal with a moderately realistic accretion 
disk—but a disk created for Interstellar by Double Negative artists rather than created 
by solving astrophysical equations such as [32] • In Appendix A.6 we give some details 


of how this and other Double Negative accretion disk images were created. This artists’ 
Interstellar disk was chosen to be very anemic compared to the disks that astronomers 



Figure 14: A moderately realistic accretion disk, created by Double Negative artists 


gravitationally lensed by the same black hole with a/M = 0.999 as in Figure 13 
with the same geometry. 


and 

and 
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see around black holes and that astrophysicists model — so the humans who travel near 
it will not get fried by X-rays and gamma-rays. It is physically thin and marginally 
optically thick and lies in the black hole’s equatorial plane. It is not currently accreting 
onto the black hole, and it has cooled to a position-independent temperature T = 4500K, 
at which it emits a black-body spectrum. 

Figure [T4| shows an image of this artists’ disk, generated with a gravitational lensing 
geometry and computational procedure identical to those for our paint-swatch disk. 
Figure 13 (no frequency shifts or associated colour and brightness changes; no lens 


flare). Christopher Nolan and Paul Franklin decided that the flattened left edge of the 
black-hole shadow, and the multiple disk images alongside that left edge, and the off- 
centred disk would be too confusing for a mass audience. So—although Interstellar^ 
black hole had to spin very fast to produce the huge time dilations seen in the movie— 
for visual purposes Nolan and Franklin slowed the spin to a/M = 0.6, resulting in the 
disk of Figure [T^. 


4-1.2. Colour and brightness changes due to freguency shifts The influences of Doppler 
and gravitational frequency shifts on the appearance of this disk are shown in Figures 
[T^,c. 

Since the left side of the disk is moving toward the camera and the right side away 
with speeds of roughly 0.55c, their light frequencies get shifted blueward on the left 
and redward on the right—by multiplicative factors of order 1.5 and 0.4 respectively 
when one combines the Doppler shift with a ~ 20 percent gravitational redshift. These 
frequency changes induce changes in the disk’s perceived colours (which we compute by 
convolving the frequency-shifted spectrum with the sensitivity curves of motion picture 


him) and also induce changes in the disk’s perceived brightness] see Appendix A.6 for 
some details. 

In Figure [T^, we have turned on the colour changes, but not the corresponding 
brightness changes. As expected, the disk has become blue on the left and red on the 
right. 

In Figure we have turned on both the colour and the brightness changes. Notice 
that the disk’s left side, moving toward the camera, has become very bright, while 
the right side, moving away, has become very dim. This is similar to astrophysically 
observed jets, emerging from distant galaxies and quasars; one jet, moving toward Earth 
is typically bright, while the other, moving away, is often too dim to be seen. 


4-2. Lens flare and the accretion disk in the movie Interstellar 

Christopher Nolan, the director and co-writer of Interstellar, and Paul Franklin, the 
visual effects supervisor, were committed to make the him as scientihcally accurate 
as possible—within constraints of not confusing his mass audience unduly and using 
images that are exciting and fresh. A fully realistic accretion disk. Figure [TSp , 
that is exceedingly lopsided, with the hole’s shadow barely discernible, was obviously 




Gravitational Lensing by Black Holes in Astrophysics and in Interstellar 


27 



14 


but with the black 


Figure 15: (a) The moderately realistic accretion disk of Figure 
hole’s spin slowed from a/M = 0.999 to a/M = 0.6 for reasons discussed in the text, (b) 
This same disk with its colours (light frequencies v) Doppler shifted and gravitationally 
shifted, (c) The same disk with its specihc intensity (brightness) also shifted in accord 
with Liouville’s theorem, ly oc. This image is what the disk would truly look like to 
an observer near the black hole. 
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Figure 16: The accretion disk of Figure 15 r (no colour or brightness shifts) with lens 
flare added—a type of lens flare called “veiling flare”, which has the look of a soft glow 
and is very characteristic of IMAX camera lenses. This is a variant of the accretion disk 
seen in Interstellar. (Figure created by our Double Negative team using DNGR, and TM 
&: ©Warner Bros. Entertainment Inc. (sl5)) This image may be used under the terms 
of the Creative Commons Attribution-NonCommercial-NoDerivs 3.0 (CC BY-NC-ND 
3.0) license. Any further distribution of these images must maintain attribution to the 
author(s) and the title of the work, journal citation and DOI. You may not use the 
images for commercial purposes and if you remix, transform or build upon the images, 
you may not distribute the modihed images. 


unacceptable. 

The hrst image in Figure [T^ the one without frequency shifts and associated colour 
and brightness changes, was particularly appealing, but it lacked one element of realism 
that few astrophysicists would ever think of (though astronomers take it into account 
when modelling their own optical instruments). Movie audiences are accustomed to 
seeing scenes hlmed through a real camera—a camera whose optics scatter and diffract 
the incoming light, producing what is called lens flare. As is conventional for movies 
(so that computer generated images will have visual continuity with images shot by 
real cameras), Nolan and Franklin asked that simulated lens flare be imposed on the 
accretion-disk image. The result, for the hrst image in Figure [Tsl is Figure 

This, with some embellishments, is the accretion disk seen around the black hole 
Cargantua in Interstellar. 

All of the black-hole and accretion-disk images in Interstellar were generated using 
DNCR, with a single exception: When Cooper (Matthew McConaughey), riding in the 
Ranger spacecraft, has plunged into the black hole Cargantua, the camera, looking 
back upward from inside the event horizon, sees the gravitationally distorted external 
universe within the accretion disk and the black-hole shadow outside it — as general 
relativity predicts. Because DNCR uses Boyer-Lindquist coordinates, which do not 
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extend smoothly through the horizon, this exceptional image had to be constructed by 
Double Negative artists manipulating DNGR images by hand. 


4-3. Some details of the DNGR accretion-disk simulations 


4-3.1. Simulating lens flare In 2002 one of our authors (James) formulated and 
perfected the following (rather obvious) method for applying lens flare to images. The 
appearance of a distant star on a camera’s focal plane is mainly determined by the point 
spread function of the camera’s optics. For Christopher Nolan’s dims we measure the 
point spread function by recording with HDR photography (see e.g. HU) a point source 
of light with the full set of 35mm and 65mm lenses typically used in his IMAX and 
anamorphic cameras, attached to a single lens reflex camera. We apply the camera’s 
lens flare to an image by convolving it with this point spread function. (For these 
optics concepts see, e.g., [M].) For the image 15i, this produces Figure [l6j More recent 
work [19] does a more thorough analysis of reflections between the optical elements in a 
lens, but requires detailed knowledge of each lens’ construction, which was not readily 
available for our Interstellar work. 


4-3.2. Modelling the accretion disk for Interstellar As discussed above, the accretion 
disk in Interstellar was an artist’s conception, informed by images that astrophysicists 
have produced, rather than computed directly from astrophysicists’ accretion-disk 
equations such as [32] . 

In our work on Interstellar, we developed three different types of disk models: 

• an inhnitely thin, planar disk, with colour and optical thickness defined by an 
artist’s image; 

• a three dimensional ‘voxel’ model; 

• for close-up shots, a disk with detailed texture added by modifying a commercial 
renderer. Mantra [50] . 

We discuss these briefly in [Appendix A.6 

5. Conclusion 

In this paper we have described the code DNGR, developed at Double Negative Ltd., for 
creating general relativistically correct images of black holes and their accretion disks. 
We have described our use of DNGR to generate the accretion-disk images seen in the 
movie Interstellar, and to explain effects that influence the disk’s appearance: light-ray 
bending, Doppler and gravitational frequency shifts, shift-induced colour and brightness 
changes, and camera lens flare. 

We have also used DNGR to explore, for a camera orbiting a fast spinning black 
hole, the gravitational lensing of a star field on the celestial sphere—including the 
lensing’s caustics and critical curves, and how they influence the stellar images’ pattern 
on the camera’s sky, the creation and annihilation of image pairs, and the image motions. 
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Elsewhere ISB we describe our use of DNGR to explore gravitational lensing by 
hypothetical wormholes; particularly, the influence of a wormhole’s length and shape on 
its lensing of stars and nebulae; and we describe the choices of length and shape that 
were made for Interstellar^ wormhole and how we generated that movie’s wormhole 
images. 
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Appendix A. Some Details of DNGR 

Appendix A. 1. Ray-tracing eguations 

In this section we give a step-by-step prescription for computing (i) the ray-tracing map 
from a point {6cs, 4>cs} on the camera’s local sky to a point {6', 0'} on the celestial sphere, 
and also (ii) the blue shift fc/ f' of light emitted from the celestial sphere with frequency 
f and received at the camera with frequency fc- Throughout we set the black hole’s 
mass M to unity, so distances are measured in terms of M. 

The foundations for our prescription are: (i) the Kerr metric written in Boyer- 
Lindquist coordinate^ 

ds^ = —a^dt^ -\- {p‘^/A)dr‘^ p^dO"^ -\- ■uj'^{d(j) — ojdtY , (A.l) 

where 

p = y/r'^ + cos^ 6 , A = — 2r , E = (r^ -|- A sin^ 6 , 

§ Our choice of Boyer-Lindquist coordinates prevents the camera from descending into the black hole. 
For that, it is not hard to switch to ingoing Kerr coordinates; but if one wants to descend even further, 
across an inner (Cauchy) horizon of the maximally extended Kerr spacetime, one must then switch 
again, e.g., to outgoing Kerr coordinates. Alain Riazuelo (private communication) has created an 
implementation of the ray tracing equations that does such multiple switches. As a foundation for 
Interstellar, Thorne did so but only for the switch to ingoing Kerr coordinates, as Interstellar assumes 
the inner horizons have been replaced by the Marolf-Ori shock singularity [S2] and the Poisson-Israel 
mass-inflation singularity [S3]; see Chapters 24 and 26 of m- 
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p-y/A 2ar S sin 9 

ce = —=— , ui = —^ , w = - 


(A.2) 


E ’ “ E2 p 

(ii) the 3+1 split of spacetime into space plus time embodied in this form of the metric; 
and (iii) the family of hducial observers (FIDOs) whose world lines are orthogonal to 
the 3-spaces of constant t, and their orthonormal basis vectors that he in those 3-spaces 
and are depicted in Figure 


0-r — 


TA d 


- 


- 

vu d(j) 


(A.3) 


p dr ’ 

We shall also need three functions of r and of a ray’s constants of motion b (axial angular 
momentum) and q (Carter constant), which appear in the ray’s evolution equations 


(A. 15) below—i.e., in the equations for a null geodesic: 


P = r‘^ + — ab , i? = — A[(6 — a)^ + g] , 0 = g — cos^ 6 


62 


sin2 0 


— a 


(A.4) 


And we shall need the function qo{bo), which identihes the constants of motion for rays 
(photons) that are unstably trapped in constant-r orbits around the black hole. This 
function is dehned parametrically by |1B] 

_r3(r3 - 6r2 + 9ro - 4a2) 

Qo — ’ ' ' ’ 9 / ’ 1 \ 9. ’ 


bo = 


' O 


3r2 + a^ro + 


a(ro-l) ’ a2(r„-l)2 

where Vo is the radius of the trapped orbit, which runs over the interval ri < ro < r 2 , 
with 


ri = 2 < 1 + cos 


2 . 
- arccos( 
o 


, r 2 = 2 1 + cos 


- arccos(+a) 
o 


.(A.6) 


The prescription that we use, in DNGR, for computing the ray-tracing map and 


the blue shift is the following concrete embodiment of the discussion in Section 2.1 


(i) Specify the camera’s location (+, 6*c, </>c), and its speed [3 and the components 
Bg, and B^ of its direction of motion relative to the FIDO at its location; and 
specify the ray’s incoming direction {9cs,4>cs) on the camera’s local sky. [Note: If 
the camera is in a circular, equatorial geodesic orbit around the black hole, then 


(A.7) 


Br = 0 , Bg = 0 , P 7 = 1 , (3 = —{fl — co) , where Q = ■ - 

« (a + rc ) 

is the geodesic angular velocity at the camera’s radius Vc, and the other quantities 


are dehned in Equations (A.2). 


(ii) Compute, in the camera’s proper reference frame, the Cartesian components (Figure 

of the unit vector N that points in the direction of the incoming ray 

= sin 6cs cos (pcs , Ny = sin 6cs sin (pcs , = cos 6cs ■ (A.8) 

(iii) Using the equations for relativistic aberration, compute the direction of motion of 
the incoming ray, np, as measured by the FIDO in Cartesian coordinates aligned 
with those of the camera: 

-VT^N, 


npy — 


-Ny + 13 

npx = 


l-l3Ny 


I- l3Ny 


npz — 


1-DN„ 


; (A.9) 
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and from these, compute the components of ni? on the FIDO’s spherical 
orthonormal basis: 


B 


npf — —- npx Bf npy -\- 


BfB^ 


K 


n 


F4i 


Bf 

= - Up X + B^ np y + 




K 

BqB^ 


K 


K 


'^F z 1 ^F§ — B§npy nnpz , 

npz , where k = Jl — B^ = jBj + i??A. 10 ) 


IV 


Compute the ray’s canonical momenta (covariant coordinate components of its 4- 
momentum) with its conserved energy —pt set to unity as a convenient convention: 

P 


Pt 1 , Pr Ep ^F f ) Pd Epp np Q , P(j) EpZU Hp ^ , 

where Ep = - (A.11) 

a + uwup^ 

is the energy measured by the FIDO. (Note: Pa can also be regarded as the ray’s 
wave vector or simply as a tangent vector to the ray.) Then compute the ray’s other 
two conserved quantities: its axial angular momentum b and its Carter constant q: 


b = p^ , q = Pg + cos 9 


sin^e 


— a 


(A. 12 ) 


V 


(Note: if we had not set pt = —1, then b would be —p^/pt and q would be (Carter 
const ant )/pj.) 

Determine, from the ray’s constants {b, q}, whether it comes from the horizon or the 
celestial sphere by the following algorithm: (a) Are both of the following conditions 
satished?: (v.l) 

bi <b <b 2 , where bi = bo{r 2 ) , 62 = bo{ri) , (A.13) 


with ri and r 2 given by Equations (A. 6 ) and the function bo{ro) given by Equations 


(A.5); and (v. 2 ) the ray’s value of q lies in the range 


q < qoib) , where qo{b) is dehned by Equations (A.5) 


(A. 14) 


(b) If the answer to (a) is yes, then there are no radial turning points for that {b, q}, 
whence if > 0 at the camera’s location, the ray comes from the horizon; and if 
Pr < 0 there, it comes from the celestial sphere, (c) If the answer to (a) is no, then 
there are two radial turning points for that {b,q}, and if the camera radius Vc is 
greater than or equal to the radius Tup of the upper turning point, then the ray 
comes from the celestial sphere; otherwise, it comes from the horizon. Here Tup is 


the largest real root of R{r) = 0, with R{r) dehned by Equations (A.4). 


VI 


If the ray comes from the celestial sphere, then compute its point of origin there 
as follows: Beginning with the computed values of the constants of motion and 
canonical momenta, and beginning in space at the camera’s location, integrate the 
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ray equation^ 
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vii 


'R + AQ' 
.. 2Ap2 ^ 


(A. 15) 

numerically, backward in time from C = 0toC = C/, where C/ is either — cx) or some 
very negative value. (This super-Hamiltonian version of the ray equations is well 
behaved and robust at turning points, by contrast, for example, with a commonly 
used version of the form p‘^dr/d( = ±a/R, p^dO/dC, = ±a/ 0 , ... ; Eqs. (33.32) of 
MTW [33].) The ray’s source point (6^', 0') on the celestial sphere is O' = 0{(f), 

0' = 0(C/)- 

Compute the light’s blue shift, i.e. the ratio of its frequency fc as seen by the camera 
to its frequency f as seen by the source at rest on the celestial sphere, from 

fc ^ (V^EE] EalE 
r 


Appendix A.2. Ray-hundle equations 


a 


(A.16) 


In this section we give a step-by-step prescription for evolving a tiny bundle of rays (a 
light beam) along a reference ray, from the camera to the celestial sphere, with the goal 
of learning the beam’s major and minor angular diameters 6^ and (5_ on the celestial 
sphere, and the angle p from e^/ on the celestial sphere to the beam’s major axis (cf. 
Figure [^. As in the previous subsection, we set the black hole’s mass to unity, M = 1. 


Our prescription is a concrete embodiment of the discussion in Section 2^ it is 
a variant of the Sachs optical scalar equations (3^; and it is a slight modification of a 
prescription developed in the 1970s by Serge Pineault and Robert Roeder 1^ 

Our prescription relies on the following functions, which are defined along the 
reference ray. (a) The components of the ray bundle’s 4-momentum (wave vector) 
on the FIDO’s orthonormal spherical basis 

^ a Pf) 1 Prh b 


p =pf = 


p 


-Pr , P = Pg 


p 


> p =p^ 


w 


W 


p^ = JpI + pI + pI ; (^•1'^) 


II These are the super-Hamiltonian variant of the null geodesic equations; see, e.g., Section 33.5 of [33] 
or Appendix A of [54| . 

^ The Pineault-Roeder equations (which were optimised for rays that begin near the black hole and 
travel from there to earth) become asymptotically singular for our situation, where the rays begin on 
the celestial sphere and travel to a camera near the black hole. We have reworked the Pineault-Roeder 
analysis so as to make their equations fully nonsingular and robust for our situation. This reworking 
entailed changing their null vector pP from 2“^/^(l, —1,0, 0) to 2“^/^(l,-|-1,0, 0) in Equation (6) of 
[35]. This simple change of one foundational equation produced many changes in the final evo lution 
equations for the ray bundle. Most importantly, it changed l/(p* + to 1/(A — p'') in Eqs. (A.19) 


and (A.21); the former diverges at the celestial sphere for our ingoing ray bundles, while the latter is 


finite there. For a detailed derivation of our prescription, see 
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M = -^ i Sa^VAr sin^ 6 cos 6 (pj. — 

4p3S2(/-pt)\ 

+2{p^ — p^) l^asin^ (pg (a — r)(a + r) cos26* + 2V^p*rsin26*^ 

+Pg — 3a^r'^ — 6r^^ — 2p^ cot 6^ j (A. 18) 

—a^p^ sin6' \^2VApq sin6' a^(r — 1) cos26' + a^{r — 1) + 2r^(r + 3) 

+a^Ap‘ cos 36* + cos 9 a^(7p* — 4p’’) + a?r{p^ (15r — 22) — Sp’’ (r — 2)) 

_4^3(pr _ 2pi (^7- _ 3))j j I , 

(c) A complex component 4/o* of the Weyl tensor, evaluated on vectors constructed from 
the ray’s 4-momentum; its real and imaginary parts work out to be: 

^ 2 (^-l)(p^-pf)2 {^^ “ 4 )' “ + ‘^P'r^^P' ^+Pi ~ 

+6p? [pj - {p^f w -p*p^S{w -1) - p^iw + 1)] + 2p^ [p^S{w - 1) (p| - 3p|) 

+3(p*)^m; + 3{p^fp^S{w - 1) - 3p‘ (w + 2) (p| -p^)] - 3(p*)^w - 2(p*)^p^ A(w - 1) 
+6(p‘)^[pi(u: + l)-p?] -2 p*p^^(m;-1)(p?-3p?)) +2 g2P0 (3p^w(p^-p|) 

+p^^(l -w) + 3pj[p^S{w -l)+p^{w + 2)] - pf [-^(w - 1) (pj - 3p|) + 3{p^fS{w - 1) 

+6p‘p^(w + 2)] + {pyS{w - 1) + 3(p*)2p^ (w + 2 ) -p*^(m; - 1 ) (pj - 3p|) ) | , (A. 19) 

- ^2 (6p|p> + P^[3p* w + p^ S{w - 1)] 

+3p? [pjiw + 1) - 2{p^fw-p^p^S{w - 1) -p?] +Pr [p^s{w - 1) (p^ - 3p|) 

+3(p*)^w + 3(p*)2p^^(M; - 1) - 3p‘ (w + 2) (p| -p^)] - (p‘)^p^^(w - 1)] 

+3(p')^ [pj-pl{w + 1)] -p^p^S{w - 1) (p| - 3p|) + Qipg [3p^w (pj -p|) (A.20) 

+p?^(l - w) + 3p^[p* A(w - 1) +p^ (w + 2)] -pf [-^(w - 1) (p| - 3p^) + 3(p*)^A(w - 1) 

+6p‘p^ {w + 2)] + {p^Siw - 1) + 3(p*)^p^ {w + 2) -p^ S{w - 1) (pj - 3p?) ) | . 

Here Qi, Q 2 , w, and S are given by 

r(r^ — 3a^ cos^ 6*) a cos 6*(3r^ — cos^ 6*) 

“ ))6 ’ ^2 — ; 

Aa^sin^6* 3ay/A{r‘^ + a!^) sinO 

(r2 + a2)2 ’ ^ S2 


w = 


(A.21) 
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We shall state our prescription for computing the shape and orientation of the 
ray bundle on the celestial sphere separately, for two cases: a ray bundle that begins 
circular at the camera; and one that begins elliptical. Then, we shall briefly sketch the 
Pineault-Roeder [35l [36] foundations for these prescriptions. 


Appendix A.2.1. Circular ray bundle at camera If the ray bundle begins at the 
camera with a tiny circular angular diameter 5cs as measured by the camera, then 
our prescription for computing its shape and orientation at the celestial sphere is this: 


(i) Introduce five new functions along the reference ray: m(C), 9{C)i h{C), and 

xiC) [which are denoted x, y, p, q and y by Pineault and Roeder]. Give them and 
their derivatives the following initial conditions at the camera, C = 0- 

M, V, g, h, V, g, h, y all vanish ; it = 1 . (A.22) 


Here a dot means d/d(. 

(ii) Integrate the following differential equations backward along the ray toward its 
source at the celestial sphere, i.e. from ( = 0 to (j. These are coupled to the ray 
equations (A. 15) for the reference ray. 


u = —'^{g cos 'ip + h sin -0) , v = —'^{g sin ip — h cos pj) , 

g = —"^{ucosp! + vsinpj) , h =— v cospj) , x = AA . (A.23) 


Here 

T = |vl/o.| = y[3fi(vl/o*)]2 + [A(vl/o*)]2, ^ = arg(To*) - 2y . (A.24) 

(iii) Upon reaching the celestial sphere at C = (f, evaluate the angular diameters of the 
bundle’s major and minor axes using the equations 

(5+ = d'{Vu‘^ + + tjg'^ + h^) , — \lg"^ + h‘^) ■ (A.25) 


Here 6' is the angular radius that the bundle would have had at the celestial sphere 
if spacetime curvature had not affected it: 

\/l — f 1 — bu^ 


5' = A 


l-/3Ny 


a 


(A.26) 


where the right hand side is to be evaluated at the camera. Also evaluate the 
bundle’s orientation angle at the celestial sphere from 


= y + ^ arg[(M + ii)){g + ih)] , 


(A.27) 


with the right hand sides evaluated at the celestial sphere, C = C/j cind with 
i = 1. The first term, y, is the contribution from parallel transporting the 

bundle’s major axis along the reference ray; the second term is the deviation from 
parallel transport caused by spacetime curvature. 
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Appendix A.2.2. Elliptical ray bundle at camera Our prescription, for a ray bundle 
that begins elliptical at the camera, is this: 

(i) In the camera’s proper reference frame, specify the Cartesian components of the 
bundle’s major axis {J^, Jy, Jz} (a unit vector), subject to the constraint that J be 


orthogonal to the ray direction n = —N [Equation (A.8)], JxUx + Jy^y + Jznz = 0. 

(ii) Compute the major-axis direction in the FIDO frame using the transformation 
equations 

■v/l /3 Jy j _ j [dnxJy J _ J (dnzJy 

^ F T. ^ T. ~ 


Jpy 


J F — Jr 


l + fdUy l + fdny l + 

where (3 is the camera’s speed (in the y direction) relative to the FIDO at its 
location. 

(iii) In the camera’s reference frame, specify (a) the bundle’s angular diameter 5cs along 
its major axis, and (b) the ratio e of the angular diameter along the minor axis to 
that along the major axis. These quantities are the same in the FIDO frame as in 
the camera frame: they are insensitive to special relativistic aberration. 

(iv) Then the initial conditions for the integration variables u, v, g, h, x are 

u,v,g,h,XiX all vanish; ft = a ■ Ji?(l-|-e)/2 , h = b • Ji?(l-|-e)/2 , 

^ = a- Ji.(l-e)/2 , h = b ■ Ji.(l - e)/2 . (A.29) 

Here a ■ = OxJfx + CLyJpy + clzJfz is the Cartesian inner product, and a and b 

are two orthogonal unit vectors that span the plane orthogonal to the hducial ray: 
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^F§^F<j> 


UFi 


bf = n 
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0 (1 


F4> 


HFi 


(v) Continue and conclude with steps (ii) and (iii) of Appendix A.2.1 


(A.30) 


Appendix A.2.3. Foundations for these ray-bundle prescriptions For completeness we 
briefly describe the Pineault-Roeder |35l [36] foundations for these ray-bundle-evolution 
prescriptions. 

At an arbitrary location C, along the reference ray, in the reference frame of a FIDO 
there, denote by {upf, the components of the hducial ray’s unit tangent 

vector. (These are the same ufA as in the previous two subsections, extended to the 
hducial ray’s entire world line.) And at this arbitrary (j, introduce the unit basis vectors 
a and b [Eqs. (A.30)] that span the plane orthogonal to the hducial ray. 

Notice that far from the black hole (at the celestial sphere), where uft = —1 
(negative because the reference ray is incoming), Upg = Up^ = 0, the only nonzero 
components of these basis vectors are = 1 and 6^ = 1; so a = and b = e^. 
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We can think of the transverse plane at ( as the complex plane, and think of the a 
direction as the real axis for complex numbers, and b as the imaginary axis. Then any 
transverse vector Y can be represented by a complex number Y, with 

Y = Y,a + nb, Y = Ya + tY,. (A.31) 

To describe the ray bundle and its evolution, (following Pineault and Roeder) we 
introduce two complex numbers (transverse vectors) ^ and p, whose real and imaginary 
parts are functions that are evolved by the ray-bundle equations (A.23): 

^ = u + iv , 7] = g + ih . (A.32) 

The outer edge of the ray bundle, at location ( along the reference ray, is described by 
the complex number 

Y = + 7]e-^^)G^ , (A.33) 

where x is also evolved by the ray-bundle equations. Here a is a parameter that ranges 
from 0 to 271, and as it changes, it carries the head of the vector Y around the bundle’s 
edge. That vector, of course, is given by 

Y = 3?[(^ + r] e-*'")e*^]a + A[(e e*'" + r] (A.34) 


[cf. Eq. (A.31)], and its tail sits on the reference ray while its head sits on the edge of 
the bundle. One can verify that the shape of this Y(ct), with a varying from 0 to 27r, is 
an ellipse. 

Let us denote by u and A the arguments of ^ and rj, so ^ = |^| rj = \ri\e^^. Then 
Eq. (A.33) reads Y = |^|e*0+°'+><:) + As a varies from 0 to 27r, the modulus 

of Y (equal to the length of the vector Y) reaches its largest value when the phase of 
the two terms is the same, i.e. when i> + a = X — a, or equivalently when a = |(A — z/). 
This maximal value of |Y| is 

|Y|max = |Y|max = |^| + Ivl , (A.35) 

and the argument of Y at the max (the angle of Y to a) is 

fi = u + a + x = ^{^ + X)+x = ^arg(^?7) + x. (A.36) 

The major diameter (not angular diameter) of the elliptical ray bundle, as measured by 
a FIDO at the location (, is this |Y|max multiplied by dp = 5cs(/c//f): 

(major diameter of ellipse) = (jF|Y|max = <^f(|^| + \v\) ■ (A.37) 

When the bundle reaches the celestial sphere, its measured angular diameter is 
the rate of increase of this major diameter with distance traveled. But at the celestial 
sphere, the distance traveled is equal to the affine parameter (, so the measured angular 
diameter is 5F(d|Y|max/dC)- Now, the FIDO-measured light frequency fp is equal to /' 
at the celestial sphere, and the real and imaginary parts of ^ = u + iv and rj = g + ih 
are increasing linearly with (, so it turns out that the formula 5F(d|Y|max/dC) becomes 
the first of Eqs. (A.25); and the angle /i of the major axis to the real axis a, Eq. ( |A.36 ), 
becomes Eq. (A.27). 
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By an argument similar to the first part of the paragraph before last, one can 
deduce that anywhere along the evolving ray bundle 

(minor diameter of ellipse) = ^^lYlmin = <5^(1^! — \ri\) , (A.38) 


and correspondingly that the minor angular diameter at the celestial sphere is the second 
of Eqs. (|A.25|). 


The evolution equations (A.23) are the equation of geodesic deviation for the 
separation vector whose spatial part in a FIDO reference frame is the vector Y [Eq. 
(A.34)] at hxed a; see Section II of Pineault and Boeder [SB], and see 


Appendix A.3. Filtering 

Appendix A.3.1. Spatial filtering and how we handle point stars In DNGR we treat 
stars as point sources of light with inhnitesimal angular size. We trace rays backwards 
from the camera to the celestial sphere. If we were to treat these rays as inhnitely thin, 
there would be zero probability of any ray intersecting a star; so instead we construct 
a beam with a narrow but hnite angular width, centred on a pixel on the camera’s sky 
and extending over a small number of adjacent pixels, we evolve the shape of this beam 
along with the ray, and if the beam intercepts the star, we collect the star’s light into 
it. This gives us some important benehts: 

• The images of our unresolved stars remain small: they don’t stretch when they get 
magnihed by gravitational lensing. 

• The fractional change in the beam’s solid angle is directly related to the optical 
magnihcation due to gravitational lensing and hence to the intensity (brightness) 
change in the image of an unresolved star. 

• When sampling images of accretion discs or extended structures such as interstellar 
dust clouds, we minimise moire artefacts by adapting the resampling hlter according 
to the shape of the beam [56] . 

Another consequence is that each star contributes intensity to several pixels. The 
eye is sensitive to the sum of these and as a single star crosses this grid, this sum can 
vary depending on the phase of the geometric image of the star on the gridj^ This can 
cause distracting flickering of stars as they traverse the virtual camera’s image plane. We 
mitigate this effect by modulating the star’s intensity across the beam with a truncated 
Gaussian hlter and by setting the beam’s initial radius to twice the pixel separation. 

With this hlter, the sum is brightest when the geometric image of the star falls 
exactly on a pixel and dimmest when centred between four pixels. The shape and width 
of the hlter are designed to give a maximum 2% diherence between these extremes, 
which we found, empirically, was not noticeable in fast changing Interstellar scenes. 


+ The same is true, of course, for each bit of a nebula or accretion disk. We discuss some details of 
how DNGR handles accretion disks in Appendix A.6 
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This result assumes the final size of the ellipse is small and the shape of the 
beam does not change significantly between adjacent pixels. In extreme cases these 
assumptions can break down, leading to a distortion in the shape of a star’s image, 
flickering, and aliasing artefacts. In these cases we can trace multiple beams per pixel 
and resample (See e.g. Chapter 7 of [S7j). However, Interstellar^ black hole, Gargantua, 
when visualized, has a moderate spin a/M = 0.6, which gives rise to few areas of extreme 
distortion, so we did not observe these problems in images created for the movie. 


Appendix A.3.2. Temporal filtering: motion blur For moving images, we need to filter 
over time as well as space. A traditional film camera typically exposes the film for half 
the time between individual frames of film, so for a typical 24 fps (frame per second) 
film, the exposure time will be 1/48 s. (This fraction is directly related to the camera’s 
“shutter angle”: the shutters in film cameras are rotating discs with a slice cut out. The 
film is exposed when the slice is over the emulsion, and the film moves onto the next 
frame when covered. The most common shape is a semicircular disk which is a “180 
degree shutter”.) Any movement of the object or camera during the exposure leads to 
motion blur. In our case, if the camera is orbiting a black hole, stellar images close 
to a critical curve appear to zip around in frantic arcs and would appear as streaks of 
light on a 1/48 s photograph. We aim to reproduce that effect when creating synthetic 
images for films. 

This motion blur in computer graphics is typically simulated with Monte Carlo 
methods, computing multiple rays per pixel over the shutter duration. In order to 
cleanly depict these streaks with Monte Carlo methods, we would need to compute the 
paths of many additional rays and the computational cost of calculating each ray is very 
high. Instead, in DNGR we take an analytic approach to motion blur (cf. Figure [Al| by 
calculating how the motion of the camera influences the motion of the deflected beam: 

The camera’s instantaneous position influences the beam’s momentary state, and 
likewise the camera’s motion affects the time derivatives of the beam’s state. We 


augment our ray and ray bundle equations (A. 15) and (A.23) to track these derivatives 


throughout the beam’s trajectory and end up with a description of an elliptical beam 
being swept through space during the exposure. The beam takes the form of a swept 
ellipse when it reaches the celestial sphere, and we integrate the contributions of all 
stars within this to create the motion-blurred imagej^ These additional calculations 
approximately double the computation time, but this is considerably faster than a naive 
Monte Carlo implementation of comparable quality. 


Appendix A.3.3. Formal description of our analytic approach to motion blur We can 
put our analytic method into a more formal setting, using the mathematical description 
of motion blurred rendering of Sung et al. [60], as follows. We introduce the quantity 


Ixy 


= V/ / r{u:H)gi{u:H)Li{u}fi)dtdiu: , 

, JnJAT 


(A.39) 


We took a similar approach when motion-blurring the accretion disks discussed in Appendix A.6 
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Figure Al: Left-to-riglit: No motion blur; Monte Carlo motion blur with four time 
samples per-pixel; Analytic motion blur. Relative computation time is approximately 
in the ratio 1:4:2. Lens flare has been added to each image to illustrate how these images 
would typically be seen in the context of a movie. 


which represents the resulting intensity of a sample located at coordinates {x, y) on the 
virtual image plane of the camera, corresponding to a ray direction u, and subtending a 
solid angle, hi. The sum is over each object, I, in the scene, AT is the time the shutter 
is open, and Li{u,t) represents the radiance of the object directed back towards the 
camera along uj from object 1. The shutter and optics may filter the incident radiance 
(specihc intensity) and this is accounted for by the term r{u,t)- The term gi{u,t) is a 
geometric term that takes the value 1 if object I is visible or 0 if occluded. 

Considering just the distortion of the celestial sphere to start with, we start with 
a ray in the direction uq at time to, in the middle of the exposure. We determine 
whether this ray gets swallowed up by the black hole or continues to the celestial sphere 
using equations (A.13) and (A.14). This corresponds to g{uo,to). We evolve the initial 
direction u and solid angle hi subtended by the sample into a new direction at the 
celestial sphere, u)', and ellipse 11', using the method described in Appendix A.l and 
Appendix A.2 


When the ray strikes the celestial sphere, we integrate over the area 11', weighting 
the integration with our truncated Gaussian hlter, corresponding to r{u). (We assume 
there is no time-dependent aspect to this hlter and ignore the t term. We also assume 
the celestial sphere is static.) 

Our motion blur method approximates L{u}) with the following expression: 


L 



(A.40) 


where the partial derivative is evaluated at the centre of the beam ooq and at to. This 
represents a ray that can be swept over an arc in the celestial sphere by varying the 
value of t. This is convolved with the beam shape 11' and hlter and occlusion terms to 
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give: 



(A.41) 


In this form, the consequences of our approximations become clearer: By evaluating 
g only at the beam centre uq and shutter-centre to, a sample that crosses the black hole’s 
shadow during the exposure would be given a constant value of g (either 0 or 1) instead 
of the correct expression that transitions during the exposure. This would lead to a 
hard edge at the shadow if the camera performed a fast pan. 

Luckily, in Interstellar^ the design of the shots meant this rarely happened: the 
main source of motion blur was the movement of the spaceship around the black hole, 
and not local camera motion. 

To handle any cases where this approximation might cause a problem, we also 
implemented a hybrid method where we launch a small number of Monte Carlo rays 
over the shutter duration and each ray sweeps a correspondingly shorter analytic path. 

We used a very similar technique for the accretion disk. 


Appendix A.4- Implementation 


DNGR was written in C-I--I- as a command-line application. It takes as input the 
camera’s position, velocity, and held of view, as well as the black hole’s location, mass 
and spin, plus details of any accretion disk, star maps and nebulae. 

Each of the 23 million pixels in an IMAX image dehnes a beam that is evolved as 
described in [Appendix A.2 


The beam’s evolution is described by the set of coupled hrst and second order 


differential equations (A.15) and (A.23) that we put into fully hrst order form and then 


numerically integrate backwards in time using a custom implementation of the Runge- 
Kutta-Fehlberg method (see, e.g.. Chapter 7 of |SH])- This method gives an estimate 
of the truncation error at every integration step so we can adapt the step size during 
integration: we take small steps when the beam is bending sharply relative to our 
coordinates, and large steps when it is bending least. We use empirically determined 
tolerances to control this behaviour. Evolving the beam along with its central ray triples 
the time per integration step, on average, compared to evolving only the central ray. 

The beam either travels near the black hole and then heads oh to the celestial 
sphere; or it goes into the black hole and is never seen again; or it strikes the accretion 
disk, in which case it gets attenuated by the disk’s optical thickness to extinction or 
continues through and beyond the disk to pick up additional data from other light 
sources, but with attenuated amplitude. We use automatic diherentiation [51] to track 
the derivatives of the camera motion through the ray equations. 

Each pixel can be calculated independently, so we run the calculations in parallel 
over multiple CPU cores and over multiple computers on our render-farm. 

We use the OpenVDB library [6T] to store and navigate volumetric data and 
Autodesk’s Maya |H2] to design the motion of the camera. (The motion is chosen to 
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fit the film’s narrative.) A custom plug-in running within Maya creates the command 
line parameters for each frame. These commands are queued up on our render-farm for 
off-line processing. 

Appendix A.5. DNGR Gode characteristics and the Double Negative render-farm 

A typical IMAX image has 23 million pixels, and for Interstellar we had to generate 
many thousand images, so DNGR had to be very efficient. It has 40,000 lines of C-I--I- 
code and runs across Double Negative’s Linux-based render-farm. Depending on the 
degree of gravitational lensing in an image jj] it typically takes from 30 minutes to several 
hours running on 10 CPU cores to create a single IMAX image. The longest renders 
were those of the close-up accretion disk when we shoe-horned DNGR into Mantra. For 
Interstellar, render times were never a serious enough issue to influence shot composition 
or action. 

Our London render-farm comprises 1633 Dell-M620 blade servers; each blade 
has two 10-core E5-2680 Intel Xeon CPUs with 156GB RAM. During production of 
Interstellar, several hundred of these were typically being used by our DNGR code. 


Appendix A. 6. DNGR modelling of accretion disks 


Our code DNGR includes three different implementations of an accretion disk: 

Thin disk: For this, we adapted our DNGR ray-bundle code so it detects 
intersections of the beam with the disk. At each intersection, we sample the artist’s 
image in the manner described in Section 2.2 and Appendix A.3 above, to determine 
the colour and intensity of the bundle’s light and we attenuate the beam by the optical 
thickness (cf. Appendix A.4). 

Volumetric disk: The volumetric accretion disk was built by an artist using SideFX 
Houdini software and stored in a volumetric data structure containing roughly 17 million 
voxels (a fairly typical number). Each voxel describes the optical density and colour of 
the disk in that region. 

We used Extinction-Based Sampling [63] to build a mipmap volume representation 
of this data. The length of the ray bundle was split into short piecewise-linear segments 
which, in turn, were traced through the voxels. We used the length of the major axis 
of the beam’s cross section to select the closest two levels in the mipmap volume; and 
in each level we sampled the volume data at steps the length of a voxel, picking up 
contributions from the colour at each voxel and attenuating the beam by the optical 
thickness. The results from the two mipmap levels were interpolated before moving on 
to the next line segment. 


ji More specifically: the render time depends mainly on the number of steps in the integration. This, 
in turn, depends on the truncation errors in the integration scheme. The tolerances are tighter on 
the position of the ray than its shape, as errors in position give rise to more noticeable artefacts. 
Empirically, renders of regions close to the black hole’s shadow are much slower than any other. 
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Close-up disk with procedural textures: Side Effects Software’s renderer, Mantra, has 
a plug-in architecture that lets you modify its operation. We embedded the DNGR ray¬ 
tracing code into a plug-in and used it to generate piecewise-linear ray segments which 
were evaluated through Mantra. This let us take advantage of Mantra’s procedural 
texture^ff] and shading language to create a model of the accretion disk with much 
more detail than was possible with the limited resolution of a voxelised representation. 
However this method was much slower so was only used when absolutely necessary. 

Disk layers close to the camera The accretion-disk images in Interstellar were 
generated in layers that were blended together to form the hnal images. Occasionally 
a layer only occupied space in the immediate vicinity of the camera, so close that 
the infuences of spacetime curvature and gravitational redshifts were negligible. These 
nearby layers were rendered as if they were in flat spacetime, to reduce computation 
time. 

Frequency shifts. In our visualizations that included Doppler and gravitational 
frequency shifts, we modelled the accretion disk as a constant-temperature blackbody 
emitter and we estimated the star temperatures from published NASA data 05 ]. 
The Doppler and gravitational frequency shifts correspond to a temperature shift in 
the blackbody spectra on arrival at the camera. We convolve those temperature- 
shifted black body spectra with the published spectral sensitivity curves of a typical 
motion picture him [HI] to generate separate red, green and blue radiance values, 
R, G and B. For the volumetric disk, these are calculated at each step and used 
in place of artist-created radiance values. In Figure [T^ we removed the doppler- 
induced intensity change by dividing Figure [TSji’s RGB triplet of radiance values 
by a weighted mean that, empirically, leaves the eye-perceived colours unchanged: 
{R',G',B'} = {R, G, R}/(0.30i?-|-0.59G-|-O.llR). This dimmed the blue side of the 
disk and brightened the red side. 

We set the white balance of our virtual camera to render a 6500 K blackbody 
spectrum with equal red, green and blue pixel values by applying a simple gain to 
each colour channel. We did not model the complex, nonlinear interaction between the 
colour-sensitive layers that occurs in real him. 


Appendix A.l. Comparison of DNGR with astrophysical codes and with film-industry 
CGI codes 


Near the end of Section ]_T, we compared our code DNGR with the state-of-the-art 
astrophysical visualization code GRay. The most important diherences—DNGR’s use 
of light-beam mappings vs. GRay’s use of individual-ray mappings, and DNGR’s use 
of ordinary processors vs. GRay’s use of GPUs—were motivated by our diherent goals: 
smoothness of images in movies vs. fast throughput in astrophysics. 


ft A procedural texture uses an algorithm, such as fractal noise, to generate detail at arbitrary resolution, 
unlike a texture based on an image file, which has a hnite resolution. 
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So far as we are aware, our code is unique in using light-beam-based mappings. 
No other code, astrophysical or him CGI, uses them. However, some him CGI codes 
use hat-spacetime mappings from the camera’s image plane to the source plane that 
are mathematically equivalent to our light beams. Specihcahy, they use ray differential 
techniques described by Igehy [50], which track, along a ray, the transverse derivatives 
of the locations of nearby rays, thereby producing the matrix \ \dff / dxf\, where 
is a ray’s source-plane location as a function of its image-plane location. Tracking 
these diherentials in hat spacetime is moderately easy because they are constant 
along straight paths between rehections, refractions, etc; and they change via simple 
matrix transformation laws at interaction locations. Tracking them through curved 
spacetime would entail the same geodesic-deviation technique as underlies our ray- 
bundle mappings. 

In Appendix A.3.1 we have described our methods of imaging star helds, using 
stars that are point sources of light: we feed each star’s light into every light beam that 
intersects the star’s location, with an appropriate weighting. 

Other gravitational-lensing codes deal with star helds diherently. For example, 
Muller and Frauendiener [25] calculate where each point-star ends up in the camera’s 
image plane, ehectively producing the inverse of the traditional ray-tracing algorithm. 
They do this for lensing by a non-spinning black hole (Schwarzchild metric) which has 
a high degree of symmetry, making this calculation tractable. Doing so for a spinning 
black hole would be much more challenging. 

A common method for rendering a star-held is to create a 2D source picture 
(environment map) that records the stars’ positions as hnite-sized dots in the source 
plane (see e.g. [HS])- This source picture is then sampled using the evolved rays from 
the camera. This has the disadvantage that stars can get stretched in an unrealistic way 
in areas of extreme magnihcation, such as near the critical curves described in Section 
3.3.1 of this paper. As we discussed in |Appendix A.3.T our DNGR light-beam technique 
with point stars circumvents this problem—as does the Muher-Frauendiener technique 
with point stars. 
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